#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
#include <time.h>

#define ITER 45  //Numero de iteracoes que serao executadas

int main (int argc, char *argv[]) 
{
	//Variaveis para medir tempo de execucao
	clock_t t0, tf;

	//Variaveis do algoritmo
	mpf_t aN, aNP1;
	mpf_t bN, bNP1;
	mpf_t tN, tNP1;
	mpf_t pN, pNP1;
	mpf_t pi;

	//Abre o arquivo onde sera escrito o resultado de pi calculado
	FILE *piFile = fopen("out_gaussLegendreSeq.txt","w+");

	if(piFile == NULL) 
	{
		return EXIT_FAILURE;
	}

	long i;

	printf("Starting...");
	printf("\n");
	//Inicializa o cronometro
	t0 = clock();
	
	//Setar precisão (usar todas as 10 milhões de casas) e inicializar variaveis do metodo
	mpf_set_default_prec(40000000);

	mpf_init(aN);
	mpf_set_ui(aN,1); //a0=1

	mpf_init(bN);
	mpf_sqrt_ui(bN,2);
	mpf_ui_div(bN,1,bN); //b0 = 1/(2^(1/2))

	mpf_init_set_str(tN,"0.25",10); //t0 = 0.25

	mpf_init(pN);
	mpf_set_ui(pN,1); //p0 = 1

	mpf_init(aNP1);
	mpf_init(bNP1);
	mpf_init(tNP1);
	mpf_init(pNP1);
	mpf_init(pi);
	
	//Loop principal
	for(i=0 ; i<ITER ; i++) 
	{   
		//Calcula an+1, bn+1, pn+1 e tn+1
		mpf_add(aNP1,aN,bN);
		mpf_div_ui(aNP1,aNP1,2); //an+1 = (an + bn)/2

		mpf_mul(bNP1,aN,bN);
		mpf_sqrt(bNP1,bNP1); // (an*bn)^(1/2)

		mpf_sub(tNP1,aN,aNP1);
		mpf_mul(tNP1,tNP1,tNP1);
		mpf_mul(tNP1,tNP1,pN);     
		mpf_sub(tNP1,tN,tNP1);  // tn+1 = tn - pn*(an-an+1)^2

		mpf_init(pNP1);
		mpf_mul_ui(pNP1,pN,2); //pn+1 = 2*pn

		mpf_set(aN,aNP1);
		mpf_set(bN,bNP1);
		mpf_set(tN,tNP1);
		mpf_set(pN,pNP1);   
	}
	
	//Depois do fim das iteracoes calcula o valor de pi
	mpf_add(pi,aN,bN);
	mpf_pow_ui(pi,pi,2);
	mpf_div_ui(pi,pi,4);
	mpf_div(pi,pi,tN); //pi = ((an + bn)^2)/(4*tn)
	
	tf = clock();
	printf("Time: %lf seconds - Iterations: %ld\n",(double)((tf - t0)/CLOCKS_PER_SEC), i);
	
	//Escreve o valor de pi no arquivo
	mpf_out_str(piFile,10,250000,pi);

	//Fecha o arquivo e limpa as variaveis do algoritmo
	fclose(piFile);
	mpf_clears(pi,aN,bN,tN,pN,aNP1,bNP1,tNP1,pNP1,NULL);
}